##################poornoobcfd - cfdnoob#####################################
###################modified from Mikael Mortensen's mesh tools code (found on github), Thank you for this!######### 
##############	mesh converter -- fluent mesh to data files of cell, face and nodes ##################



import re
import sys
from sets import Set
from copy import copy
from numpy import zeros, array,savetxt
import scitools
import scitools.filetable

try: 
  infile = sys.argv[1]
except:
  print 'usage : ', sys.argv[0], 'fluent_mesh.msh'



re_dimline   = re.compile(r"\(2\s(\d)\)")
re_comment   = re.compile(r"\(0\s.*")
re_zone_init = re.compile(r"\(10\s\(0\s(\w+)\s(\w+)\s(\d+)\s(\d+)\)\)")
re_zone      = re.compile(r"\(10\s\((\w+)\s(\w+)\s(\w+)\s(\d+)\s(\d)\)(\(|)")
re_face_init = re.compile(r"\(13(\s*)\(0\s+(\w+)\s+(\w+)\s+(0|0 0)\)\)")
re_face      = re.compile(r"\(13(\s*)\((\w+)\s+(\w+)\s+(\w+)\s+(\w+)\s+(\w+)\)(\s*)(\(|)")
re_periodic  = re.compile(r"\(18.*\((\w+)\s+(\w+)\s+(\w+)\s+(\w+)\).*\(")
re_pfaces    = re.compile(r"((^\s)|)(\w+)(\s*)(\w+)")
re_cells_init= re.compile(r"\(12(\s*)\(0(\s+)(\w+)(\s+)(\w+)(\s+)(0|0 0)\)\)")
re_cells     = re.compile(r"\(12.*\((\w+)\s+(\w+)\s+(\w+)\s+(\d+)\s+(\d+)\)\)")
re_cells2    = re.compile(r"\(12(\s*)\((\w+)\s+(\w+)\s+(\w+)\s+(\w+)\s+(\w+)\)(\s*)(\(|)")
re_zones     = re.compile(r"\((45|39)\s+\((\d+)\s+(\S+)\s+(\S+).*\)\((.*|[0-9]+[\.]*[0-9]*)\)\)")
re_parthesis = re.compile(r"(^\s*\)(\s*)|^\s*\)\)(\s*)|^\s*\(\s*)")
re_parant	 = re.compile(r"(^\s*\)(\s*)|^\s*\)\)(\s*)|^\s*\(\s*)")


nodenos=[];xnode=[];ynode=[];znode=[];

cell_map = {}               # Maps cell id with nodes
cell_face_map = {}          # Maps cell id with faces
face_cell_map = {}          # Maps face id with two cells
face_list = []              # List of faces [[id, 2-4 nodes, 2 connecting cells and type]]
face_map = {}               # For each cell a dictionary with key=face and val=local face number
nodes = None                # Will be numpy array of nodes
xnode = []
ynode = []
znode = []
nodenos =[]

# Information about connectivity and boundaries
boundary_map = {}           # For each cell and local face number, get connected cell and local face number
boundary_val_map = {}       # Holds boundary conditions, like V, v, P, etc
boundary_nodes = {}         # List of nodes attached to a boundary. Key is zone id
boundary_nodes_face_map = {}# Map of the faces that belongs to a node on a boundary
boundary_faces = {}            # List of cells attached to a boundary facet. Key is zone id
boundary_cells = {}
# Some global values
num_cells = {}              # Total number of cells in different zones
zones = {}                  # zone information
zone_number_of_faces = {}   # number of faces for each zone

def read_periodic(ifile, periodic_dx):
        """Scan past periodic section. Periodicity is computed by FEniCS."""
        while 1:
            line = ifile.readline()
            a = re.search(re_pfaces, line)
            if a:
                continue
            break

def read_zone_vertices(dim, Nmin, Nmax, ifile):
        """Scan ifile for vertices and add to mesh_editor."""
        # First line could be either just "(" or a regular vertex.
        # Check for initial paranthesis. If paranthesis then read a new line, else reset
        pos = ifile.tell()
        line = ifile.readline()
        if not re.search(re_parthesis, line): 
            ifile.seek(pos) # reset            
        # read Nmax-Nmin vertices
        for i in range(Nmin, Nmax + 1):
            line = ifile.readline()
            vertex = [eval(x) for x in line.split()]
            if dim == 2:
                #nodes.append(i - Nmin, vertex[0], vertex[1])
                nodenos.append(i-Nmin),xnode.append(float(vertex[0])),ynode.append(float(vertex[1]))
            else:
                #nodes.append(i - Nmin, vertex[0], vertex[1], vertex[2])
                nodenos.append(i-Nmin),xnode.append(float(vertex[0])),ynode.append(float(vertex[1])),
                znode.append(float(vertex[2]))

def read_faces1(zone_id, Nmin, Nmax, bc_type, face, ifile):
        """Read all faces and create cell_map + boundary maps."""        
        pos = ifile.tell() # current position
        line = ifile.readline()
        if not re.search(re_parthesis, line): # check for initial paranthesis. If paranthesis then read a new line, else reset
            ifile.seek(pos)
            
        # read Nmax-Nmin faces
        for i in range(Nmin, Nmax + 1):
            line = ifile.readline()
            ln = line.split()
            if face == 0:
                nd = int(ln[0], 16) # Number of vertices
                nds = [int(x, 16) for x in ln[1:(nd + 1)]]
                cells = [int(x, 16) for x in ln[(nd + 1):]]
            else:
                nd = face
                nds = [int(x, 16) for x in ln[:nd]]
                cells = [int(x, 16) for x in ln[nd:]]
                        
            if min(cells) == 0: # A boundary zone
                if zone_id in boundary_cells:
                    boundary_cells[zone_id][max(cells)] = array(nds)
                else:
                    boundary_cells[zone_id] = {max(cells): array(nds)}

            for c in cells:
                if c > 0:    
                    if not c in cell_map:
                        cell_map[c] = copy(nds)                    
                    else:
                        cell_map[c] = list(Set(cell_map[c] + nds))


def read_faces(zone_id, Nmin, Nmax, bc_type, face, lines):
    """Read all faces and create cell_face_map + some boundary maps."""
    
    if re.search(re_parant, lines[0]): # check for initial paranthesis
        dummy = lines.pop(0)
    ls = []
    for i in range(Nmin, Nmax + 1):
        line = lines[i - Nmin]
        ln = line.split()
        if face == 0:
            nd = int(ln[0]) # Number of nodes
            nds = [int(x, 16) for x in ln[1:(nd + 1)]]
            cells = [int(x, 16) for x in ln[(nd + 1):]]
        else:
            nd = face
            nds = [int(x, 16) for x in ln[:nd]]
            cells = [int(x, 16) for x in ln[nd:]]
            
        face_list.append([nd, copy(nds), copy(cells), bc_type, zone_id])
        if len(nds) == 2:
            face_cell_map[(nds[0], nds[1])] = copy(cells)
            face_cell_map[(nds[1], nds[0])] = copy(cells)

        face_number = len(face_list)
        if min(cells) == 0: # A boundary zone
            if zone_id in boundary_nodes:
                boundary_nodes[zone_id] += nds
                boundary_faces[zone_id] += [face_number - 1]
                for nd in nds:
                    if nd in boundary_nodes_face_map[zone_id]:
                        boundary_nodes_face_map[zone_id][nd] += [face_number - 1]
                    else:
                        boundary_nodes_face_map[zone_id][nd] = [face_number - 1]
            else:
                boundary_nodes[zone_id] = nds
                boundary_faces[zone_id] = [face_number - 1]
                boundary_nodes_face_map[zone_id] = { nd: [face_number - 1]}

        for c in cells:
            if c > 0:                                                                      
                if not c in cell_face_map:
                    cell_face_map[c] = [face_number]
                else:
                    # Preliminary cell_face_map. Needs shaping up later
                    cell_face_map[c].append(face_number)
                    
    if min(cells) == 0:
        boundary_nodes[zone_id] = list(Set(boundary_nodes[zone_id]))

def create_cell_face_map(dim, mesh_format):
    """Creates maps from cells to faces and gives faces local numbering."""
    if mesh_format == 'fenics':
        create_cell_map(dim)
    else:
        eval('create_cell_face_map_' + str(dim) + 'D()') 


def scan_fluent_mesh1(ifile):  
        """Scan fluent mesh and generate maps."""
        dim = 0
        one = 0
        while 1:
            line = ifile.readline()
            if len(line) == 0:
                print 'Finished reading file\n'
                break

            if dim == 0: # Dimension usually comes first
                a = re.search(re_dimline, line)
                if a: 
                    print 'Reading dimensions\n'
                    dim = int(a.group(1))
                    #editor.open(mesh, dim, dim)
                    continue
            
            if one == 0: # The total number of vertices
                a = re.search(re_zone_init, line)
                if a:
                    print 'Reading zone info\n'
                    one, num_vertices, dummy1, dummy2 = int(a.group(1)), \
                        int(a.group(2), 16), int(a.group(3), 16), int(a.group(4))
                    #editor.init_vertices(num_vertices)
                    continue
                
            a = re.search(re_zone, line) # Vertices
            if a:
                zone_id, first_id, last_id = int(a.group(1), 16), \
                    int(a.group(2), 16), int(a.group(3), 16)
                print 'Reading ', last_id - first_id + 1,' vertices in zone ', zone_id + 1, '\n'
                read_zone_vertices(dim, first_id, last_id, ifile)
                continue
                
            a = re.search(re_zones,line) # Zone info
            if a:
                print 'Reading zone info ', line
                dummy, zone_id, zone_type, zone_name, radius =  \
                        int(a.group(1)), int(a.group(2)),  a.group(3), \
                        a.group(4), a.group(5)
                zones[zone_id] = [zone_type, zone_name, radius]
                continue
            
            a = re.search(re_cells_init, line) # Get total number of cells/elements
            if a:
                print 'Reading cell info ', line
                first_id, tot_num_cells = int(a.group(3),16), int(a.group(5), 16)
                #editor.init_cells(tot_num_cells)
                continue

            a = re.search(re_cells,line) # Get the cell info.
            if a:
                zone_id, first_id, last_id, bc_type, element_type = \
                    int(a.group(1),16), int(a.group(2), 16), int(a.group(3), 16), \
                    int(a.group(4), 16), int(a.group(5), 16)
                print 'Found ', last_id - first_id + 1,' cells in zone ', zone_id, '\n'
                if last_id == 0:
                    raise TypeError("Zero elements!")
                continue

            a = re.search(re_cells2,line) # Get the cell info.
            if a:
                raise TypeError("Wrong cell type. Can only handle one single cell type")

            a = re.search(re_face_init, line)
            if a:
                print 'Reading total number of faces\n', line
                continue
                
            a = re.search(re_face, line)
            if a:
                print 'Reading faces ', line
                zone_id, first_id, last_id, bc_type, face_type = \
                    int(a.group(2), 16), int(a.group(3), 16), int(a.group(4), 16), \
                    int(a.group(5), 16), int(a.group(6), 16)
                read_faces1(zone_id, first_id, last_id, bc_type, face_type, ifile)
                continue
            
            a = re.search(re_periodic, line)
            if a:
                print 'Scanning past periodic connectivity\n', line
                read_periodic(ifile, periodic_dx)
                continue
            
            if any([re.search(st, line) for st in (re_parthesis, re_comment)]) or \
                                                                not line.strip():
                continue
                    
            # Should not make it here
            raise IOError('Something went wrong reading fluent mesh.')


def scan_fluent_mesh(lines):
    """Scan fluent mesh and generate numerous maps."""
    # Warning! Not yet tested for multiple interior zones
    dim = 0
    one = 0
    num_faces = 0
    while 1:
        try:
            line = lines.pop(0)
        except:
            print 'Finished reading file\n'
            break
        if dim == 0: # Dimension usually comes first
            a = re.search(re_dimline, line)
            if a: 
                print 'Reading dimensions\n'
                dim = int(a.group(1))
                print 'Mesh is ' + str(dim) + 'D\n'
                continue
        
        if one == 0: # The total number of nodes
            a = re.search(re_zone_init, line)
            if a:
                print 'Reading zone info\n'
                one, num_vertices, dummy1, dummy2 = int(a.group(1)), \
                     int(a.group(2), 16), int(a.group(3), 16), int(a.group(4))
                continue
            
        # a = re.search(re_zone, line) # Nodes
        # if a:
        #     zone_id, first_id, last_id, dummy1, dummy2 = int(a.group(1), 16), \
        #         int(a.group(2), 16), int(a.group(3), 16), int(a.group(4)), \
        #         int(a.group(5))
        #     print 'Reading nodes in zone ', zone_id + 1, '\n'
        #     # read_zone_vertices(dim, first_id, last_id, ifile)
        #     read_zone_nodes(dim, first_id, last_id, lines)
        #     lines = lines[(last_id - first_id + 1):]  
        #     continue
            
        a = re.search(re_zones,line) # Zone info
        if a:
            print 'Reading zone ', line
            dummy, zone_id, zone_type, zone_name, radius =  \
                       int(a.group(1)), int(a.group(2)),  a.group(3), \
                       a.group(4), a.group(5)
            zones[zone_id] = [zone_type, zone_name, radius]
            continue
        
        a = re.search(re_cells_init, line) # Get total number of cells/elements
        if a:
            print 'Reading cell info ', line
            first_id, tot_num_cells = int(a.group(3),16), int(a.group(5), 16)
            continue

        a = re.search(re_cells,line) # Get the cell info.
        if a:
            zone_id, first_id, last_id, bc_type, element_type = \
                int(a.group(1),16), int(a.group(2), 16), int(a.group(3), 16), \
                int(a.group(4), 16), int(a.group(5), 16)
            print 'Reading cells in zone ', zone_id, '\n'
            if last_id == 0:
                raise TypeError("Zero elements!")
            num_cells[zone_id] = [first_id, last_id, bc_type, element_type]
            continue

        a = re.search(re_cells2,line) # Get the cell info.
        if a:
            raise TypeError("Wrong cell type. Can only handle one single cell type")

        a = re.search(re_face_init, line)
        if a:
            print 'Reading total number of faces\n', line
            num_faces = int(a.group(3),16)
            continue
        
        
        a = re.search(re_face, line)
        if a:
            print 'Reading faces ', line
            zone_id, first_id, last_id, bc_type, face_type = \
                 int(a.group(2), 16), int(a.group(3),16), int(a.group(4),16), \
                 int(a.group(5), 16), int(a.group(6), 16)
            read_faces(zone_id, first_id, last_id, bc_type, face_type, lines)
            lines = lines[(last_id - first_id + 1):]
            zone_number_of_faces[zone_id] = last_id - first_id + 1




flumesh =open(infile,'r')
scan_fluent_mesh1(flumesh)

flumesh.close()


flumesh = open(infile,'r')
lines = flumesh.readlines()
scan_fluent_mesh(lines)
flumesh.close()

g = open('cell2node.dat','w')

tot_num_cells = len(cell_map)

for i in range(tot_num_cells):
	g.write(repr(i+1))
	for n in cell_map[i+1]:
		g.write(' %s'%(repr(n)))
	g.write('\n')

g.close()

g = open('node.dat','w')
tot_nodes = len(xnode)
nodenos=array(nodenos)
xnode=array(xnode);ynode=array(ynode);znode=array(znode)
scitools.filetable.write_columns(g,nodenos,xnode,ynode,znode)
g.close()

lfmap = len(cell_face_map)
g = open('cell2face.dat','w')
for i in range(lfmap):
	g.write(repr(i+1))
 	for n in cell_face_map[i+1]:
 		g.write(' %s'%(repr(n)))
 	g.write('\n')
g.close()


aa = zeros((len(face_list),8))



for i in range (len(face_list)):
    aa[i][0] = i+1
    aa[i][1] = face_list[i][1][0]
    aa[i][2] = face_list[i][1][1]
    aa[i][3] = face_list[i][1][2]
    aa[i][4] = face_list[i][1][3]
    aa[i][5] = face_list[i][2][0]
    aa[i][6] = face_list[i][2][1]
    aa[i][7] = face_list[i][4]

# b = aa.tolist() 

# g  = open('face2node.dat','w')
savetxt('face2node.dat',aa,fmt='%.1d',delimiter=' ',newline='\n')




